Destruction of Superconductivity by Impurities in the Attractive Hubbard Model 
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We study the effect oi U — impurities on the superconducting and thermodynamic properties of 
the attractive Hubbard model on a square lattice. Removal of the interaction on a critical fraction 
of /crit ~ 0.30 of the sites results in the destruction of off-diagonal long range order in the ground 
state. This critical fraction is roughly independent of filling in the range 0.75 < p < 1.00, although 
our data suggest that fait might be somewhat larger below half-filling than at p = 1. We also find 
that the two peak structure in the specific heat is present at / both below and above the value 
which destroys long range pairing order, ft is expected that the high T peak associated with local 
pair formation should be robust, but apparently local pairing fluctuations are sufficient to generate 
a low temperature peak. 

PACS numbers: 



The study of the interplay of disorder and interac- 
tions is one of the central problems in condensed matter 
physics, and has been analyzed through a wide range of 
techniques, d H H S H lE Quantum Monte Carlo 
(QMC) is one approach which has been applied only 
relatively recently, but now a number of investigations 
of the disordered, repulsive Hubbard model exist, both 
on finite lattices |3, la Q| and within dynamical mean 
field theory where the momentum dependence of the self- 
energy is neglected. ^3 Several interesting effects, includ- 
ing the possibility of the enhancement of TNeei by site 
disorder [ij, the occurrence of a Mott transition away 
from half-filling, "12] and a change in the temperature de- 
pendence of the resistivity from insulating to metallic be- 
havior when interactions are turned on,[l3 have been ob- 
served. These studies have also explored different types 
of disorder, including randomness in the chemical poten- 
tial, hopping, Zeeman fields, and on-site interaction. [l^ 
At half-filling, some of these disorder types preserve 
particle-hole symmetry which has an important effect 
both on the behavior of the static (magnetic and charge) 
correlations,!^ and on the conductivity. 

The work described above concerns the interplay of 
randomness and repulsive interactions. Superconductor- 
insulator transitions provide a motivation to explore 
models with an effective attractive electron-electron inter- 
action, since these novel transitions are widely studied ex- 
perimentally, especially in two dimensions . [la IT^ ITM Il9l | 
However, QMC studies of models like the attractive Hub- 
bard Hamiltonian with disorder are less numerous than 
for the repulsive case.|0|23 Previous work focused on 
locating the critical site disorder strength for the insulat- 
ing transition and determining the value of the conductiv- 
ity and its possible universalitv. [2 L] and on the interplay 
between the loss of phase coherence and the closing of 
the gap in the density of states. p3.l23l | 

Another particularly interesting application of the at- 



tractive Hubbard model with disorder is the question 
of how inhomogeneities and pairing affect each other in 
high-temperature superconductors. The attractive Hub- 
bard model is the simplest Hamiltonian where this inter- 
play can be studied qualitatively. 

Rather than looking at random chemical potentials, as 
in the previous work just described, here we will instead 
explore the effect of disordered interaction, specifically 
turning off the attraction on a fraction / of the sites in 
the —U Hubbard model. i24| Such vacancies are analo- 
gous to non-magnetic impurities in the repulsive Hub- 
bard model, and, indeed, at half-filling, the two prob- 
lems are related by a particle-hole transformation. We 
determine the filling dependence of the critical impurity 
concentration for the destruction of superconductivity, 
through measurements of the equal time pair correla- 
tions and the current-current correlations. We also study 
the effects of the disappearance of superconductivity on 
the specific heat. Without disorder, a high temperature 
peak, associated with pair formation and a low temper- 
ature one, associated with pair coherence, are present. 
With disorder, we show that the high temperature peak 
is robust as superconductivity disappears, and the low 
temperature one is weakened. 



I. THE MODEL AND COMPUTATIONAL 
APPROACH 

The attractive Hubbard Hamiltonian we study is, 

(rr')(T r<T 

- ^[/(r)(nrT-i)(nrx-i)- 

r 

Here cl^{cj.a) are fermion creation(destruction) opera- 
tors at site r with spin cr, and Ura- = cl^Cra- We use 
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the coefficient t of the kinetic energy term to set our en- 
ergy scale (i.e. t = 1). The kinetic energy lattice sum 
(rr') is over nearest neighbor sites on a two dimensional 
square lattice, and fi is the chemical potential. The on- 
site attraction U{r) is chosen to take on the two values 
U{r) =0 and t/(r) = — \U\ with probabilities / and 1 — / 
respectively. In this paper, we focus on a single in- 
teraction strength \U\ = 4. Notice that we have chosen 
a particle-hole symmetric form of the interaction. As a 
consequence, the energy level of a singly occupied va- 
cancy site is higher by \U\/2 than a singly occupied site 
with nonzero interaction. We will report on results for 
the non-symmetric form elsewhere. 

We begin by briefly reviewing the physics of the pure 
attractive Hubbard model, / = 0. At half-filling, charge 
density wave and superconducting correlations are de- 
generate, and this symmetry of the order parameters 
drives the superconducting (and CDW) transition tem- 
peratures to zero. Doping breaks the symmetry in favor 
of pairing correlations, and there is now a Kosterlitz- 
Thouless transition to a superconducting state at finite 
temperature. It is well established that Tc rises rapidly 
with doping away from half-filling, reaching a maximal 
value a,t p = 1±S with S 1/8. Tc decreases gradually 
thereafter. |23| There is some debate __as to the exact 
value of the maximal Tc, with estimates 28] in the range 
0.05^ < Tcmax < 0.2t. The problem lies with the rather 
large system sizes needed to study Kosterlitz-Thouless 
transitions numerically, and hence to benchmark the ac- 
curacy of the approximate analytic calculations. Our fo- 
cus here will not be on this finite temperature phase tran- 
sition, but rather on the quantum phase transition out of 
the superconducting ground state which occurs through 
the introduction of a non-zero fraction / of vacancies. 

To distinguish a superconducting phase, we first look 
for long range structure in the equal time pair-pair cor- 
relation function, 



c(r) = 
At(r) = 



;A(r)At(0); 



This pair-pair correlation function can also be summed 
spatially to define the structure factor. 



P,=^c(r)=^(A(r)At(0); 



The finite size scaling of Pg is described in the following 
section. 

The current-current correlations probe the superfiuid 
weight and provide an alternate means to detect the de- 
struction of superconductivity. We define, 

A,,(r,T) = (j,(r,T)j,(0,0)) 



jxirr) 



Mr 



and the Fourier transform in space and imaginary time, 
P^U^.^n) = V /''dre^'»-'-e-'""^A,,(r,r), 

r -^0 

where LOn — inir/ (3. 

The longitudinal part of the current-current correla- 
tion function satisfies the f-sum rule, which relates its 
value to the kinetic energy Kx, 

= lim^^_o Axx{qx,qy = 0,w„ = 0) 
A^ = Kx. 

Here Kx = {-tJ2ai4+i^aCv,a + claCr+x^a))^ Meanwhile, 
in the superconducting state the transverse part. 



A^ 



linig^^o Axxiqx = 0,qy,UJn = 0) 



can differ from the longitudinal part, the difference being 
the superfiuid stiffness Ds, 

Ds/n - [A^-A^] 
= [Kx - A^] . 

Thus the current-current correlations provide an alterna- 
tive, complementary method to the equal time pair cor- 
relations for looking at the superconducting transition. 

In order to evaluate these quantities, we use the de- 
terminant QMC method. [3Qi] In this approach, we dis- 
cretize the inverse temperature /3 — LAt, and a path 
integral expression is written down for the partition func- 
tion Z. The electron-electron interactions are decoupled 
by the introduction of a Hubbard-Stratonovich field. The 
fermion degrees of freedom can then be integrated out 
analytically, leaving an expression for Z which involves 
an integral over the Hubbard-Stratonovich field, with an 
integrand which is the product of two determinants of 
matrices of dimension the system size. We perform the 
integral stochastically. In the case of the attractive Hub- 
bard model considered here, the traces over the spin up 
and spin down electrons are given by the determinant of 
the same matrix, the integrand is a perfect square, and 
hence there is no sign problem. 

Observables are evaluated as the appropriate combina- 
tions of Greens functions which are given by matrix ele- 
ments of the inverse of the matrix appearing as the inte- 
grand. Systematic errors in the pairing and current mea- 
surements associated with the discretization of /3, with 
our choice of At, are typically smaller than the error 
bars associated with the statistical fiuctuations for a sin- 
gle disorder realization and the error bars associated with 
sample-to-sample variations. The 'Trotter errors' are dis- 
cernable in the energy, but we have verified they do not 
affect the qualitative behavior of the specific heat. Our 
procedure for measuring the specific heat is described in 
a subsequent section. 

As remarked above, there is no 'sign problem' for this 
attractive Hubbard model, so we can do computations 
at arbitrarily low temperature (large /3). In practice, we 
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FIG. 1: The equal time pair correlation c(r) is shown as a 
function of the separation |r| of the points of injection and 
removal of the pair. Here the lattice size is N—8x8, f — 1/16, 
and the densities p — 0.750 (left) and p — 0.875 (right). Long 
range correlations build up as the temperature is lowered, 
even though the attractive interactions have been somewhat 
diluted. Error bars are the size of the symbols. Scatter in 
the data is associated with the anisotropy in the correlations, 
which are not a function of the distance only. 



find that for the parameter values and lattice sizes un- 
der consideration here the superconducting correlations 
reach their asymptotic {T — 0) values when T < 1/12 
(/3 > 12), though we collect data also at T = 1/16 to 
check this conclusion. These temperatures are about one 
hundredth of the bandwidth W = 8t. 



II. EQUAL TIME PAIRING CORRELATIONS 

We first show the spatial behavior of the pairing cor- 
relations c(r). In Fig. ^ the vacancy fraction / = 1/16 
is fixed, and we study how long range correlations in c(r) 
develop as the temperature T is reduced. Results for two 
different densities p = 0.750 and p = 0.875 are given. The 
lattice size iV=8x8. In Fig. El we fix T = 1/16 and ex- 
amine c(r) for different /. As the vacancy rate increases, 
the long range order vanishes. 

These results can be analyzed more carefully by per- 
forming a finite size scaling study of the pair structure 
factor Pg. If the pair correlations c(r) extend over the 
entire lattice, Pg, which sums this quantity, will grow 
linearly with size. On the other hand, if the pair correla- 
tions are finite in range, Ps is independent of size. Fig.|21 
shows Ps as a function of /3 for different lattice sizes. In 
Fig. 13 we choose a small vacancy fraction / = 1/16 (left 
panels) and see that at large f3 the pair structure fac- 
tor increase significantly with lattice size. In the right 



FIG. 2: The equal time pair correlation c(r) is shown as a 
function of the separation |r| of the points of injection and re- 
moval of the pair. Here the lattice size is A'^=8x8, the density 
p — 0.750 (left), and the temperature T = 1/16. Long range 
correlations are decreased with increasing vacancy fraction /. 
Density p — 0.875 is at right. 

panels of Fig. O at large vacancy fractions, / = 4/16 
and / = 5/16, the pair structure factor grows much less 
rapidly with lattice size at low temperatures. 

Fig. ^determines the critical vacancy fraction through 
a finite size scaling analysis. As shown by Huse,[33| the 
spin-wave correction to the pair structure factor is ex- 
pected to be proportional to the linear lattice size: 



where Aq is the superconducting gap function at zero 
temperature, and a{U, f) is independent of L. In Fig. ^ 
we plot the low temperature value of Pg versus 
for lattice sizes ranging from Af=8x8 to A^= 14x14, and 
p = 0.750, p = 0.875, and p = 1.000. At smah /, there 
is the expected nearly linear behavior, extrapolating to 
a nonzero value in the thermodynamic limit. [3l| As / 
increases, the nonzero extrapolation disappears. This 
figure contains one of the central results of this paper, 
namely the determination of the critical values of va- 
cancy fraction / for the destruction of superconductivity. 
To within the resolution of these simulations, the critical 
dilution fraction /dit ~ 0.30 for the destruction of super- 
conductivity in the ground state is the same for fillings 
p = 0.750 and p ^ 0.875. 

Half-filling appears to behave somewhat differently. 
Disorder first enhances superconductivity, as seen by the 
crossing of the / = and / = 1/16 lines in Fig. 4, be- 
fore driving it to zero a bit sooner than for the doped 
cases. Impurity sites break the special superconducting 
and charge density wave degeneracy at p = 1, thereby 
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FIG. 3: (a) The equal time pair structure factor Ps is shown 
as a function of inverse temperature Here the density p = 
0.750, and dilution fraction / = 1/16. Pa grows significantly 
with lattice size A'^ at low temperatures. The temperature at 
which data for the difi'erent lattice sizes separate indicates the 
point at which the coherence length ^{T) becomes comparable 
to the linear lattice size, (b) p = 0.875, / — 1/16. (c) p — 
0.750, / = 4/16. (d) p = 0.875, / = 5/16. ' 
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FIG. 4: These finite size scaling plots show the equal time 
pair structure factor Ps at low temperature T — 1/16 as a 
function of linear lattice size 1/L = 1/%/iV. Different curves 
are for different values of the vacancy fraction. For small /, Ps 
extrapolates to a nonzero value in the thermodynamic limit 
1/L 0. For larger / there is no longer a nonzero extrapo- 
lated value. The curves are upper left panel p = 0.750, lower 
left panel p = 0.875, and upper right panel p = 1.000. An 
interesting feature of the half-filling data is the corssing of the 
f — and / = 1/16 curves, indicating an initial enhancement 
of pairing by the impurity sites. 



providing a possible reason for the initial enhancement. 

It should be noted that the values of /crit indicate that 
the transition is neither of the classical site-percolation 
type, for which jciassicai,sito ^ q nor it is related to 
the quantum percolation one, for which jquantum,sitc _ 
jquantum.bond ^ ^^le latter reflecting the fact that any 
disorder localizes tight-binding electrons, [s^ 



III. CURRENT CORRELATIONS AND 
SUPERFLUID WEIGHT 

These conclusions are confirmed by studying the 
current-current correlations. We first show, in the left 
panels of Fig. [SJ that the longitudinal current-current 
correlation Axxiqx,qy = 0,iw„ — 0) extrapolates to the 
kinetic energy as 0. This sum rule is satisfied in all 
parameter regimes, both small and large vacancy rate /, 



and both small and large temperatures T. 

In contrast, the transverse response Axx{Qx = 
Q,qy,iijJn = 0), given in the right hand panels of Fig. [SJ 
extrapolates to the kinetic energy with qy only when the 
temperature is high at small /. As pairing correlations 
develop across the lattice with decreasing T (Fig. ^ the 
transverse response breaks away from the kinetic energy, 
indicating a nonzero superfluid weight Dg. At vacancy 
concentrations beyond the point at which superconduc- 
tivity is destroyed, remains zero with decreasing tem- 
perature. 

In Fig. we show the finite-size extrapolated values 
of Ps inferred from Fig. ^ These measurements give an 
indication of the location of /crit in the thermodynamic 
limit. To within our numerical uncertainties, /crit is the 
same for the fillings p — 0.750 and p — 0.875, and may be 
a bit less for p = 1.000. In Fig. El for density p = 0.875, 
we also show the low temperature values of the superfluid 
fraction obtained from Fig.|Sl The two measurements. 
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FIG. 5: Left panels: Longitudinal current-current correlations 
^xx{<lx,qy = 0, iu>„ = 0) as a function of q^- Here the density 
p = 0.875. Axx{qx,qy = 0,iwn = 0) extrapolates to the ki- 
netic energy (symbols at g = 0) for all temperatures and both 
fillings when / = 0.04 (top), which is superconducting, and 
/ — 0.32 (bottom), which is not. Right panels: Transverse 
current-current correlations Axxiqx ~ 0,qy,iu!„ = 0) extrap- 
olate to the kinetic energy (solid square at qx — 0) at high 
temperature for / — 0.04 (top), but break away at low T 
indicating the presence of a non-zero superfluid density. For 
/ = 0.32 (bottom) the dilution fraction exceeds fait so now 
at low temperature there is no superfluid density. All data 
are for 10x10 lattices. 



Ao and Dg, give consistent results for /crit- 



IV. SPECIFIC HEAT 

At half-filling, the specific heat of the repulsive (at- 
tractive) Hubbard Hamiltonians, in the absence of impu- 
rities, exhibit two features. There is a peak 
associated with moment (pair) formation at high tem- 
peratures, T « U/3. At lower temperatures, T « 
J/4 = t'^/U, there is a second peak associated with 
magnetic ordering (pair coherence). Interestingly, in 
two dimensions, this two peak structure appears to sur- 
vive even down to weak coupling where the two en- 
ergy scales merge, l23l in contrast to the behavior in one 
dimension, ISTl l38l Is^ l40| and in the paramagnetic 
phase in infinite dimensions. |4ll l43l | 

What happens to this behavior when disorder is intro- 
duced, specifically when U = sites are inserted, as in 



1.2 
0.8 

0.6 
0.4 
0.2 

0, 



I 


•-• [ A (f) / A (f=0) f 


p=0.750 




H[A„(f)/A„(f=0) f 


p=0.875 




♦♦[A„(f)/A„(f=0) f 


p=1.000 




▲ ▲ D^(f)/D(f=0) 


p=0.875 


1 


4- iW 





0.1 



0.2 
f 



0.3 



0.4 



FIG. 6: Values for the L —> 00 pair structure factor versus 
vacancy fraction / at density p — 0.750 (circles), p = 0.875 
(squares) and p = 1.000 (diamonds). The inverse tempera- 
ture /3 = 10. Also shown is the superfluid fraction Da on a 
fixed 10x10 lattice size. Both measurements give a consistent 
estimate of fc ~ 0.30, 



our present model? To address this question we compute 
the specific heat by fitting QMC data for E^mciTn) to 
the functional form. 



M 



-/3/A 



The parameters A and c/ are selected to minimize 

{Efit{Tn) — EqniciTn))'^ 



X 



1 



Nt 

E 

n=l 



(<5i^qmc(T„))2 



We choose a number of parameters M equal to about 
one-fourth of the QMC data points to allow a good fit to 
the data without overfitting. We check that different M 
around this value all provide comparable results. 

Fig.|7|summarizes our results for the specific heat. At a 
vacancy fraction / = 0.40, the high temperature peak in 
the specific heat, which is associated with the formation 
of local pairs, is, as expected, somewhat broadened rela- 
tive to the clean system (/ = 0.00), but remains other- 
wise robust. The low temperature peak, associated with 
pair coherence, is much more substantially reduced, and 
shifted to lower T . However, it appears to remain present 
even for this value of / which is beyond /crit- 



V. CONCLUSIONS 

In this paper we have looked at the destruction of the 
superconducting state in the attractive Hubbard model 
through the systematic inclusion oiU — Q sites. A con- 
sistent picture emerges from an analysis of the pairing 
structure factor and the superfiuid fraction, namely that 
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FIG. 7: Specific heat of tfie attractive Hubbard model on 
iV^lOxlO lattices and density p = 1.000 for / = 0.00 (solid 
curve) and / = 0.40 (dashed curve). The high temperature 
peak indicates the temperature of local pair formation and 
is, as expected, not very much affected by the dilution of 
the attractive interaction on some of the sites. The low T 
peak associated with pair coherence is substantially reduced 
by dilution, but appears to be present even at large /, where 
long range order has been lost. 



pairing survives out to a dilution of about 3/10 of the at- 
tractive sites, and that this value is not very dependent 
on the filling, although it may be somewhat increased 
away from half-filling. As the fraction / of J7 = sites is 
increased the high temperature peak in C'(T), which sig- 
nals local pair formation, remains relatively unchanged. 
The low T peak which signals the establishment of pair- 
ing order gets pushed down. Although the inclusion of 
non-interacting sites might remind one of percolation, the 
transition here is not percolative. The critical value of / 
is neither that of classical site-percolation nor of quan- 
tum percolation. Finally, we have shown that directly 
at half-filling, pairing correlations are enhanced by the 
inclusion of vacancy sites. We attribute this effect to the 
breaking of degeneracy of the CDW and superconducting 
order in favor of pairing. 
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